## World War 2 Appendix Dataset ##
## Prepared by Kevin Fahey ##
## And Doug Atkinson ##

##### Clear Everything #####

rm(list=ls())

## Set Seed ##
set.seed(20120630)

##### Load in Packages #####

library(foreign)
library(sandwich)
library(stargazer)
library(tidyverse)
library(haven)
library(plm)
library(pglm)
library(lmtest)
library(stringr)
library(ggplot2)
library(maptools)
library(Matching)
library(foreach)
library(fixest)
options(scipen=7)
options(digits=4)

##### Load in ClusterMod for Std. Errors #####
clusterMod<-function(model, cluster)
{
  require(multiwayvcov)
  require(lmtest)
  vcovCL<-cluster.vcov(model, cluster)
  
  coef<-coeftest(model, vcovCL)
  #w<-waldtest(model, vcov = vcovCL, test = "F")
  get_confint<-function(model, vcovCL){
    t<-qt(.975, model$df.residual)
    ct<-coeftest(model, vcovCL)
    cse<-sqrt(diag(vcovCL))
    est<-cbind(ct[,1], cse,ct[,1]-t*ct[,2], ct[,1]+t*ct[,2])
    colnames(est)<-c("Estimate", "Clustered SE","LowerCI","UpperCI")
    return(est)
  }
  ci<-round(get_confint(model, vcovCL),4)
  return(list(coef, ci))
}

## Set Working Directory here! ##
## Read in Datasets ##

dat <- read.csv("dat.csv")
ww2.soldiers <- read.csv("ww2soldiers.csv")
final.statedat <- read.csv("finalstatedat.csv")
vitals <- read.csv("mortality rates.csv")

# Create dataset for analysis
dat.analysis <- dat[which(dat$year >= 41),]

## Create annual datasets ##
dat42 <- dat.analysis %>% filter(year == 42)
dat43 <- dat.analysis %>% filter(year == 43)
dat44 <- dat.analysis %>% filter(year == 44)
dat45 <- dat.analysis %>% filter(year == 45)

## Create matched dataset ##
## Load in Census Data ## 
## Note: to obtain this dataset, you need to obtain ICPSR 02896/DS0032/02896-0032-Data.dta
## This is publicly available on the ICPSR dataverse website
## URL: https://www.icpsr.umich.edu/web/ICPSR/studies/2896

dat.census <- merge(dat, census1940.1, by.x = "fips", by.y = "fips", all.x = T)


####################
##### ANALYSIS #####
####################

##### FIGURES #####
## Figure 4
hist(dat.analysis$propelig, col = "gray", breaks = 60,
     main = "", xlab = "Distribution of county-year proportion of enlisted individuals, World War 2.")
dev.copy(png, 'histdvappendix.png')
dev.off()


## Figure 5
hist(dat.analysis$DemAvg, col = "lightblue", breaks = 60,
     main = "", xlab = "Distribution of county-year Democratic vote share, averaged from 1932, 1936, and
1940 elections.")
dev.copy(png, 'histivappendix.png')
dev.off()


## Figure 6
hist(dat.analysis$demcong1940, col = "blue2", breaks = 60,
     main = "", xlab = "Average Democratic vote share per county-year, 1932, 1926, and 1940 Presidential elections")
dev.copy(png, 'histcongappendix.png')
dev.off()


## Figure 7 -- draftees versus volunteers
draft.durat <- ww2.soldiers[which(ww2.soldiers$drafted == 1),] # 6190144 obs
vol.durat <- ww2.soldiers[which(ww2.soldiers$volunteer == 1),] # 1750566 obs

draft.out <- draft.durat %>%
  group_by(year, month) %>%
  summarise(n()) # 375 obs
vol.out <- vol.durat %>%
  group_by(year, month) %>%
  summarise(n()) # 312 obs

draft.out <- draft.out[which(draft.out$month >= 1 & draft.out$month <= 12),] #  72 obs
vol.out <- vol.out[which(vol.out$month >= 1 & vol.out$month <= 12),] # 72 obs
draft.out$date <- (draft.out$year * 100) + draft.out$month
vol.out$date <- (vol.out$year * 100) + vol.out$month

draft.out <- draft.out[,c("n()", "date")]
vol.out <- vol.out[,c("n()", "date")]

colnames(draft.out) <- c("totdraft", "date")
colnames(vol.out) <- c("totvol", "date")

time.out <- merge(draft.out, vol.out, by.x = "date",
                  by.y = "date")
time.out$time <- seq(from = 1, to = 72)

# And plot
plot(time.out$time, time.out$totdraft, type = "l", col = "chocolate",
     lwd = 3, xlab = "Time (in Months) from January 1940",
     ylab = "Number of Enlistments", main = "",
     xaxt = "s")
lines(time.out$totvol, type = "l", col = "purple", lwd = 3, lty = 2)
legend(47, 270000, legend = c("Drafted", "Volunteers"),
       col = c("chocolate", "purple"),
       lty = 1:2, cex = 1.1)
dev.copy(png, 'draftvolrates.png')
dev.off()


## Figure 8
# Identify vars with imbalance on which to match
t.test.out <- matrix(NA, nrow = 400, ncol = 6)

for(i in 1:400){
  
  tmp <- as.data.frame(cbind(dat.census$propelig, 
                             dat.census$swingcounty,
                             dat.census[,i+85]))
  if(is.numeric(tmp$V3) == F) next
  if(min(na.omit(tmp$V3)) == max(na.omit(tmp$v3))) next
  
  test <- t.test(tmp$V3 ~ tmp$V2)
  
  t.test.out[i,2] <- as.numeric(paste(test$estimate[1], sep = ""))
  t.test.out[i,3] <- as.numeric(paste(test$estimate[2], sep = ""))
  t.test.out[i,4] <- as.numeric(paste(test$statistic, sep = ""))
  t.test.out[i,5] <- as.numeric(paste(test$p.value, sep = ""))
  t.test.out[i,6] <- as.numeric(paste(test$parameter, sep = ""))
  
  cat("\r", i, "of", 400)
}

t.test.out[,1] <- colnames(dat.census)[86:485]
t.test.out <- as.data.frame(t.test.out)
colnames(t.test.out) <- c("test.var", "group1", "group2", "tstat", "pvalue", "df")
t.test.out$pvalue <- as.numeric(paste(t.test.out$pvalue, sep = ""))

## Find where pvalue <= 0.05 ##
different.counties <- t.test.out[which(t.test.out$pvalue <= 0.05),]

## create dataset of theoretically-derived vars from different.counties ##
prematch.dat <- dat.census[,c("propelig", "swingcounty", "electionyear",
                              "GOPAvg", "pcturban", "totpop.thou",
                              "avgcropvalue", "pctmale.elig2", "pctblack",
                              "farmprop", "year", "tsch1820", "tsch2124",
                              "m25edu0", "m25edu14", "m25col13", "m25col4",
                              "rurfarm", "pctm14lf", "mfarm", "mlabor",
                              "wholwage.y", "servwage.y", "pcturb40",
                              "medperre", "pctradio", "pctrefri", "avrent",
                              "meschm25", "unempnem", "density", "fips")]

# Define vars
Y = prematch.dat$propelig
Tr = prematch.dat$swingcounty
X = prematch.dat[,c(3:32)]

t <- na.omit(cbind(Y, Tr, X))

## Estimate balance on Tr ##
tr.predict <- glm(Tr ~ ., data = prematch.dat[, -c(1, 32)])
summary(tr.predict)

Z <- cbind(t[,c(4:6, 8:11, 16:24, 26, 27, 30)])

## Match ##
match <- Match(Y=t$Y, Tr=t$Tr, X=t[,c(3:31)], 
               Z=Z, BiasAdjust=T, estimand="ATT", 
               M=1, replace=T)

# check for balance 
postmatch.dat <- as.data.frame(rbind(t[match$index.treated,],
                                     t[match$index.control,]))
tr.postdict <- glm(Tr ~ ., data = postmatch.dat[, -c(1)])

## Create covariate balance plot
library(cobalt); library(MatchIt)
balance.table <- bal.tab(match, treat = t$Tr, covs = t[,c(3:31)],
                         estimand = "ATT")
covs <- t[,c(3:31)]

# and plot
m.out <- matchit(f.build("Tr", covs), data = t,
                 method = "nearest", replace = T)
love.out <- love.plot(m.out, thresholds = c(m = 0.1, ks = 0.1),
                      stats = c("m", "ks.statistics"), binary = "std", abs = T,
                      var.order = "unadjusted", sample.names = c("Unmatched", "Matched"),
                      position = "top", shapes = c("circle filled", "diamond filled"), 
                      colors = c("mediumorchid4", "goldenrod2"))
png("loveplot.png", res = 100)
print(love.out)
dev.off() # not working; manually exported as .png


## Figure 9: Placebo Test
iter.placebo <- matrix(NA, nrow = 100, ncol = 3)
iter.placebo.42 <- matrix(NA, nrow = 100, ncol = 3)
graph.placebo <- NA
graph.placebo.42 <- NA

for(i in 10:90){
  
  tmp <- dat.analysis[,c("propelig", "swingcounty", "DemAvg", "GOPAvg",
                         "pcturban", "totpop.thou", "avgcropvalue",
                         "pctmale.elig2", "pctblack", "farmprop", "year")]
  tmp$placebo <- ifelse(tmp$DemAvg >= i & tmp$DemAvg < (i+1), 1, 0)
  tmp42 <- tmp[which(tmp$year == 42),]
  tmp42$placebo <- ifelse(tmp42$DemAvg >= i & tmp42$DemAvg < (i+1), 1, 0)
  
  model <- plm(propelig ~ placebo
               + GOPAvg 
               + pcturban + totpop.thou
               + avgcropvalue + pctmale.elig2
               + pctblack + farmprop
               + as.factor(year), data = tmp,
               effect = "time", index = c("year"),
               model = c("random"))
  model42 <- lm(propelig ~ placebo
                + GOPAvg 
                + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, data = tmp42)
  
  iter.placebo[i,1] <- model$coefficients[2]
  iter.placebo[i,2] <- sqrt(diag(vcovHC(model)))[2]
  iter.placebo[i,3] <- model$coefficients[2]/sqrt(diag(vcovHC(model)))[2]
  
  iter.placebo.42[i,1] <- model42$coefficients[2]
  iter.placebo.42[i,2] <- sqrt(diag(vcovHC(model42)))[2]
  iter.placebo.42[i,3] <- model42$coefficients[2]/sqrt(diag(vcovHC(model42)))[2]
  
  graph.placebo[i] <- sum(na.omit(tmp$placebo))
  graph.placebo.42[i] <- sum(na.omit(tmp42$placebo))
  
  cat(i, "\r")
}

# Name names #
iter.placebo <- as.data.frame(iter.placebo)
colnames(iter.placebo) <- c("coef", "se", "tstat")
iter.placebo.42 <- as.data.frame(iter.placebo.42)
colnames(iter.placebo.42) <- c("coef", "se", "tstat")

## And identify the number of stat-sig iterations ##
iter.placebo$statsig <- ifelse(abs(iter.placebo$tstat) >= 1.96, 1, 0)
sum(na.omit(iter.placebo$statsig))
iter.placebo$statsig <- as.factor(iter.placebo$statsig)
iter.placebo.42$statsig <- ifelse(abs(iter.placebo.42$tstat) >= 1.96, 1, 0)
sum(na.omit(iter.placebo.42$statsig))
iter.placebo.42$statsig <- as.factor(iter.placebo.42$statsig)

## Plot the distribution of placebo tests ##
plot(iter.placebo$tstat, pch = c(20, 17)[iter.placebo$statsig],
     main = "", ylab = "Estimated T-Statistics", xlab = "Binned Mean Vote Share for Roosevelt, 1932-1940 Elections",
     ylim = c(-5,5), xlim = c(20, 96), col = c("blue", "red")[iter.placebo$statsig])
abline(a = 1.96, b = 0, col = "black", lwd = 4, lty = 2)
abline(a = -1.96, b = 0, col = "black", lwd = 4, lty = 2)
abline(v = 50, col = "darkgrey", lwd = 4, lty = 3)
dev.copy(png, 'placeboappendix1.png')
dev.off()

## Plot the distribution of placebo tests ##
plot(iter.placebo.42$tstat, pch = c(20, 17)[iter.placebo.42$statsig],
     main = "", ylab = "Estimated T-Statistics", xlab = "Binned Mean Vote Share for Roosevelt, 1932-1940 Elections, 1942 Only",
     ylim = c(-5,5), xlim = c(20, 96), col = c("blue", "red")[iter.placebo.42$statsig])
abline(a = 1.96, b = 0, col = "black", lwd = 4, lty = 2)
abline(a = -1.96, b = 0, col = "black", lwd = 4, lty = 2)
abline(v = 50, col = "darkgrey", lwd = 4, lty = 3)
dev.copy(png, 'placeboappendix2.png')
dev.off()

## Coefficient Plot of Placebos ##
library(arm)
coefs <- iter.placebo.42$coef[28:90]
sds <- iter.placebo.42$se[28:90]
lo.ci <- coefs - (sds * 1.96)
hi.ci <- coefs + (sds * 1.96)
names <- seq(from = 28, to = 90, by = 1)

plot(coefs ~ names, pch = 20, col = "darkgrey", lwd = 3, main = "",
     ylab = "Coefficients, Placebo Test Regression Estimates",
     xlab = "Bin Starting Location (Bin = N + 1)")



## Figure 10 -- Randomization Inference
# Create Output Matrix
n <- 2000
random.inference <- matrix(NA, ncol = 3, nrow = n)

for(i in 1:n){
  
  tmp <- dat[,c("fips", "year", "propelig", "swingcounty", 
                "GOPAvg", "pctmale.elig2", "totpop.thou",
                "pcturban", "pctblack", "avgwage",
                "avgcropvalue", "farmprop")]
  tmp$placebo <- sample(tmp$swingcounty, size = length(tmp$swingcounty), replace = F)
  
  model <- plm(propelig ~ placebo + GOPAvg 
               + pcturban + pctmale.elig2
               + totpop.thou + avgcropvalue 
               + pctblack + farmprop
               + as.factor(year), data = tmp, 
               effect = "time", index = c("year"),
               model = c("random"))
  
  random.inference[i,1] <- model$coefficients[2]
  random.inference[i,2] <- sqrt(diag(vcovHC(model)))[2]
  random.inference[i,3] <- model$coefficients[2]/sqrt(diag(vcovHC(model)))[2]
  
  cat(i, "\r") 
}

random.inference <- as.data.frame(random.inference)
colnames(random.inference) <- c("coef", "se", "tstat")

# And how many are statistically-significant
random.inference$statsig <- ifelse(abs(random.inference$tstat) >= 1.96, 1, 0)
sum(random.inference$statsig)
random.inference$statsig <- as.factor(random.inference$statsig)

# How many in expected direction
neg <- sum(ifelse(random.inference$tstat <= -1.96, 1, 0))
post <- sum(ifelse(random.inference$tstat >= 1.96, 1, 0))

# Plot the distribution of placebo tests
hist(random.inference$tstat, col = "darkgrey", breaks = 80,
     main = "", xlab = "t-statistics (Swing county coefficient/standard error)")
abline(v = -1.96, lwd = 4, col = "black")
abline(v = 1.96, lwd = 4, col = "black")
dev.copy(png, filename = 'randominfhist.png')
dev.off()



##### TABLES #####
## Table 1 -- Correlation plot, DVs ##
cor_dvs <- dat.analysis %>% dplyr::select(propelig, propdraft,
                                   male1554.elig, draft_melig)
xtable(cor(cor_dvs, use = "pairwise.complete.obs"), digits = 2)


## Manuscript Table 2 -- Main Models ##

# original/baseline models
main_all <- feols(propelig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop
                  + yr42 + yr43 + yr44 + yr45,
                  data = dat.analysis,
                  cluster = "fips")
main_elec <- feols(propelig ~ swingcounty + electionyear
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop,
                   data = dat.analysis,
                   cluster = "fips")
main_inter <- feols(propelig ~ swingcounty * electionyear
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
main_42 <- feols(propelig ~ swingcounty
                 + GOPAvg + pcturban + totpop.thou
                 + avgcropvalue + pctmale.elig2
                 + pctblack + farmprop, 
                 data = dat42,
                 se = "hetero")
main_43 <- feols(propelig ~ swingcounty
                 + GOPAvg + pcturban + totpop.thou
                 + avgcropvalue + pctmale.elig2
                 + pctblack + farmprop, 
                 data = dat43,
                 se = "hetero")
main_44 <- feols(propelig ~ swingcounty
                 + GOPAvg + pcturban + totpop.thou
                 + avgcropvalue + pctmale.elig2
                 + pctblack + farmprop, 
                 data = dat44,
                 se = "hetero")
main_45 <- feols(propelig ~ swingcounty
                 + GOPAvg + pcturban + totpop.thou
                 + avgcropvalue + pctmale.elig2
                 + pctblack + farmprop, 
                 data = dat45,
                 se = "hetero")
main_am1 <- feols(propelig ~ splitswing
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop
                  + yr42 + yr43 + yr44 + yr45, 
                  data = dat.analysis,
                  cluster = "fips")
main_am2 <- feols(propelig ~ pickupswing
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop
                  + yr42 + yr43 + yr44 + yr45, 
                  data = dat.analysis,
                  cluster = "fips")

# draftees only
draft_all <- feols(propdraft ~ swingcounty
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45,
                   data = dat.analysis,
                   cluster = "fips")
draft_elec <- feols(propdraft ~ swingcounty + electionyear
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
draft_inter <- feols(propdraft ~ swingcounty * electionyear
                     + GOPAvg + pcturban + totpop.thou
                     + avgcropvalue + pctmale.elig2
                     + pctblack + farmprop,
                     data = dat.analysis,
                     cluster = "fips")
draft_42 <- feols(propdraft ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat42,
                  se = "hetero")
draft_43 <- feols(propdraft ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat43,
                  se = "hetero")
draft_44 <- feols(propdraft ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat44,
                  se = "hetero")
draft_45 <- feols(propdraft ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat45,
                  se = "hetero")
draft_am1 <- feols(propdraft ~ splitswing
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")
draft_am2 <- feols(propdraft ~ pickupswing
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")

# male-eligible interpolated DV
melig_all <- feols(male1554.elig ~ swingcounty
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45,
                   data = dat.analysis,
                   cluster = "fips")
melig_elec <- feols(male1554.elig ~ swingcounty + electionyear
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
melig_inter <- feols(male1554.elig ~ swingcounty * electionyear
                     + GOPAvg + pcturban + totpop.thou
                     + avgcropvalue + pctmale.elig2
                     + pctblack + farmprop,
                     data = dat.analysis,
                     cluster = "fips")
melig_42 <- feols(male1554.elig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat42,
                  se = "hetero")
melig_43 <- feols(male1554.elig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat43,
                  se = "hetero")
melig_44 <- feols(male1554.elig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat44,
                  se = "hetero")
melig_45 <- feols(male1554.elig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat45,
                  se = "hetero")
melig_am1 <- feols(male1554.elig ~ splitswing
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")
melig_am2 <- feols(male1554.elig ~ pickupswing
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")

# male-eligible interpolated DV, draftees only
delig_all <- feols(draft_melig ~ swingcounty
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45,
                   data = dat.analysis,
                   cluster = "fips")
delig_elec <- feols(draft_melig ~ swingcounty + electionyear
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop,
                    data = dat.analysis,
                    cluster = "fips")
delig_inter <- feols(draft_melig ~ swingcounty * electionyear
                     + GOPAvg + pcturban + totpop.thou
                     + avgcropvalue + pctmale.elig2
                     + pctblack + farmprop,
                     data = dat.analysis,
                     cluster = "fips")
delig_42 <- feols(draft_melig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat42,
                  se = "hetero")
delig_43 <- feols(draft_melig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat43,
                  se = "hetero")
delig_44 <- feols(draft_melig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat44,
                  se = "hetero")
delig_45 <- feols(draft_melig ~ swingcounty
                  + GOPAvg + pcturban + totpop.thou
                  + avgcropvalue + pctmale.elig2
                  + pctblack + farmprop, 
                  data = dat45,
                  se = "hetero")
delig_am1 <- feols(draft_melig ~ splitswing
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")
delig_am2 <- feols(draft_melig ~ pickupswing
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + yr42 + yr43 + yr44 + yr45, 
                   data = dat.analysis,
                   cluster = "fips")


# E-Tables
etable(main_all, main_elec, main_inter, 
       main_42, main_43, main_44, main_45, 
       main_am1, main_am2,
       cluster = "fips", tex = TRUE,
       drop = "fips", digits = "r3", digits.stats = "r3",
       title = "Impact of county partisanship on U.S. Army enlistment rates in World War 2. Unit-clustered standard errors listed in parentheses.",
       label = "reg1",
       fitstat = c("n", "ar2", "f"),
       dict = c(swingcounty = "Swing County (50-55)",
                GOPAvg = "GOP Vote Share",
                pcturban = "Pct Urban",
                totpop.thou = "Cty Pop (1000s)",
                avgcropvalue = "Avg Crop Value",
                pctmale.elig2 = "Pct Male Eligible",
                pctblack = "Pct Black",
                farmprop = "Prop. Farms",
                yr42 = "1942", yr43 = "1943", yr44 = "1944", yr45 = "1945",
                electionyear = "Election Year",
                splitswing = "Swing County (47.5-52.5)",
                pickupswing = "Swing County (45-50)"))

etable(draft_all, draft_elec, draft_inter, 
       draft_42, draft_43, draft_44, draft_45, 
       draft_am1, draft_am2,
       cluster = "fips", tex = TRUE,
       drop = "fips", digits = "r3", digits.stats = "r3",
       title = "Impact of alternate specifications of swing-county status on U.S. Army enlistment rates in World War 2. Year and state fixed effects employed in all models. Unit- clustered fixed effects employed in all models.",
       label = "robust1",
       fitstat = c("n", "ar2", "f"),
       dict = c(swingcounty = "Swing County (50-55)",
                GOPAvg = "GOP Vote Share",
                pcturban = "Pct Urban",
                totpop.thou = "Cty Pop (1000s)",
                avgcropvalue = "Avg Crop Value",
                pctmale.elig2 = "Pct Male Eligible",
                pctblack = "Pct Black",
                farmprop = "Prop. Farms",
                yr42 = "1942", yr43 = "1943", yr44 = "1944", yr45 = "1945",
                electionyear = "Election Year",
                splitswing = "Swing County (47.5-52.5)",
                pickupswing = "Swing County (45-50)"))

etable(melig_all, melig_elec, melig_inter, 
       melig_42, melig_43, melig_44, melig_45, 
       melig_am1, melig_am2,
       cluster = "fips", tex = TRUE,
       drop = "fips", digits = "r3", digits.stats = "r3",
       title = "Impact of alternate specifications of swing-county status on U.S. Army enlistment rates in World War 2. Year and state fixed effects employed in all models. Unit- clustered fixed effects employed in all models.",
       label = "robust1",
       fitstat = c("n", "ar2", "f"),
       dict = c(swingcounty = "Swing County (50-55)",
                GOPAvg = "GOP Vote Share",
                pcturban = "Pct Urban",
                totpop.thou = "Cty Pop (1000s)",
                avgcropvalue = "Avg Crop Value",
                pctmale.elig2 = "Pct Male Eligible",
                pctblack = "Pct Black",
                farmprop = "Prop. Farms",
                yr42 = "1942", yr43 = "1943", yr44 = "1944", yr45 = "1945",
                electionyear = "Election Year",
                splitswing = "Swing County (47.5-52.5)",
                pickupswing = "Swing County (45-50)"))

etable(delig_all, delig_elec, delig_inter, 
       delig_42, delig_43, delig_44, delig_45, 
       delig_am1, delig_am2,
       cluster = "fips", tex = TRUE,
       drop = "fips", digits = "r3", digits.stats = "r3",
       title = "Impact of alternate specifications of swing-county status on U.S. Army enlistment rates in World War 2. Year and state fixed effects employed in all models. Unit- clustered fixed effects employed in all models.",
       label = "robust1",
       fitstat = c("n", "ar2", "f"),
       dict = c(swingcounty = "Swing County (50-55)",
                GOPAvg = "GOP Vote Share",
                pcturban = "Pct Urban",
                totpop.thou = "Cty Pop (1000s)",
                avgcropvalue = "Avg Crop Value",
                pctmale.elig2 = "Pct Male Eligible",
                pctblack = "Pct Black",
                farmprop = "Prop. Farms",
                yr42 = "1942", yr43 = "1943", yr44 = "1944", yr45 = "1945",
                electionyear = "Election Year",
                splitswing = "Swing County (47.5-52.5)",
                pickupswing = "Swing County (45-50)"))


## Manuscript Figure 5 -- Aggregation to State ##
statem1 <- feols(enlistpop ~ swing32
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem2 <- feols(enlistpop ~ swing36
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem3 <- feols(enlistpop ~ swing40
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem4 <- feols(enlistpop ~ swingavg
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")

statem5 <- feols(enlistmale ~ swing32
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem6 <- feols(enlistmale ~ swing36
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem7 <- feols(enlistmale ~ swing40
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem8 <- feols(enlistmale ~ swingavg
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")

statem9 <- feols(draftpop ~ swing32
                 + as.factor(year)
                 + as.factor(statefips), 
                 data = final.statedat, se = "hetero")
statem10 <- feols(draftpop ~ swing36
                  + as.factor(year)
                  + as.factor(statefips), 
                  data = final.statedat, se = "hetero")
statem11 <- feols(draftpop ~ swing40
                  + as.factor(year)
                  + as.factor(statefips), 
                  data = final.statedat, se = "hetero")
statem12 <- feols(draftpop ~ swingavg
                  + as.factor(year)
                  + as.factor(statefips), 
                  data = final.statedat, se = "hetero")

statem13 <- feols(draftmale ~ swing32
                  + as.factor(year)
                  + as.factor(statefips), 
                  data = final.statedat, se = "hetero")
statem14 <- feols(draftmale ~ swing36
                  + as.factor(year)
                  + as.factor(statefips), 
                  data = final.statedat, se = "hetero")
statem15 <- feols(draftmale ~ swing40
                  + as.factor(year)
                  + as.factor(statefips), 
                  data = final.statedat, se = "hetero")
statem16 <- feols(draftmale ~ swingavg
                  + as.factor(year)
                  + as.factor(statefips), 
                  data = final.statedat, se = "hetero")

etable(statem1, statem2, statem3, statem4, statem5, statem6, statem7, statem8,
       cluster = "statefips", tex = TRUE,
       drop = c("statefips", "year"), digits = "r3", digits.stats = "r3",
       title = "Impact of state presidential vote share on U.S. Army enlistment rates in World War 2. Year and state fixed effects employed in all models. Unit- clustered standard employed in all models.",
       label = "stateagg",
       fitstat = c("n", "ar2", "f"),
       dict = c(swing32 = "1932 Election",
                swing36 = "1936 Election",
                swing40 = "1940 Election",
                swingavg = "Avg. of Elections")
       )

etable(statem9, statem10, statem11, statem12, statem13, statem14, statem15, statem16,
       cluster = "statefips", tex = TRUE,
       drop = c("statefips", "year"), digits = "r3", digits.stats = "r3",
       title = "Impact of state presidential vote share on U.S. Army enlistment rates in World War 2. Year and state fixed effects employed in all models. Unit- clustered standard employed in all models.",
       label = "stateagg",
       fitstat = c("n", "ar2", "f"),
       dict = c(swing32 = "1932 Election",
                swing36 = "1936 Election",
                swing40 = "1940 Election",
                swingavg = "Avg. of Elections")
)

## Congressional Results -- Figure 6 in Manuscript ##
mcavg <- lm(propelig ~ swingcongavg
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat)
se.mcavg <- sqrt(diag(vcovHC(mcavg)))

mc36a <- lm(propelig ~ swingcong36
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat42)
se.mc36a <- sqrt(diag(vcovHC(mc36a)))

mc36b <- lm(propelig ~ splitcongswing36
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat42)
se.mc36b <- sqrt(diag(vcovHC(mc36b)))

mc36c <- lm(propelig ~ pickupcongswing36
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat42)
se.mc36c <- sqrt(diag(vcovHC(mc36c)))

mc38a <- lm(propelig ~ swingcong38
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat42)
se.mc38a <- sqrt(diag(vcovHC(mc38a)))

mc38b <- lm(propelig ~ splitcongswing38
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat42)
se.mc38b <- sqrt(diag(vcovHC(mc38b)))

mc38c <- lm(propelig ~ pickupcongswing38
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat42)
se.mc38c <- sqrt(diag(vcovHC(mc38c)))

mc40a <- lm(propelig ~ swingcong40
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat42)
se.mc40a <- sqrt(diag(vcovHC(mc40a)))

mc40b <- lm(propelig ~ splitcongswing40
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat42)
se.mc40b <- sqrt(diag(vcovHC(mc40b)))

mc40c <- lm(propelig ~ pickupcongswing40
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat42)
se.mc40c <- sqrt(diag(vcovHC(mc40c)))

mc42a <- lm(propelig ~ swingcong42
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat43)
se.mc42a <- sqrt(diag(vcovHC(mc42a)))

mc42b <- lm(propelig ~ splitcongswing42
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat43)
se.mc42b <- sqrt(diag(vcovHC(mc42b)))

mc42c <- lm(propelig ~ pickupcongswing42
            + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat43)
se.mc42c <- sqrt(diag(vcovHC(mc42c)))

# And output
library(stargazer)
stargazer(mcavg, mc36a, mc36b, mc36c, mc38a, mc38b, mc38c,
          se = list(se.mcavg, se.mc36a, se.mc36b, 
                    se.mc36c, se.mc38a, se.mc38b, se.mc38c),
          no.space = T, df = F,
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          title = "Effect of Congressional swing-county status on enlistment rates in World War 2",
          dep.var.labels = c("Enlisted (100s)"))

stargazer(mc40a, mc40b, mc40c, mc42a, mc42b, mc42c,
          se = list(se.mc40a, se.mc40b, se.mc40c, 
                    se.mc42a, se.mc42b, se.mc42c),
          no.space = T, df = F,
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          title = "Effect of Congressional swing-county status on enlistment rates in World War 2",
          dep.var.labels = c("Enlisted (100s)"))


## Regional Analysis -- Table 2 in Manuscript ## 
mreg1 <- plm(propelig ~ swingcounty + newengland + confederacy
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop
             + as.factor(year),
             data = dat[which(dat$year >= 41),],
             effect = "time", index = c("year"),
             model = c("random"))
se.mreg1 <- coeftest(mreg1, 
                     vcov=vcovHC(mreg1, type = "sss", cluster = "time"))[,2]

mreg2 <- plm(propelig ~ swingcounty
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop
             + as.factor(region)
             + as.factor(year),
             data = dat[which(dat$year >= 41),],
             effect = "time", index = c("year"),
             model = c("random"))
se.mreg2 <- coeftest(mreg2, 
                     vcov=vcovHC(mreg2, type = "sss", cluster = "time"))[,2]

mreg42 <- lm(propelig ~ swingcounty + newengland + confederacy
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop,
             data = dat[which(dat$year == 42),])
se.mreg42 <- sqrt(diag(vcovHC(mreg42)))

mreg43 <- lm(propelig ~ swingcounty + newengland + confederacy
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop,
             data = dat[which(dat$year == 43),])
se.mreg43 <- sqrt(diag(vcovHC(mreg43)))

mreg44 <- lm(propelig ~ swingcounty + newengland + confederacy
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop,
             data = dat[which(dat$year == 44),])
se.mreg44 <- sqrt(diag(vcovHC(mreg44)))

mreg45 <- lm(propelig ~ swingcounty + newengland + confederacy
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop,
             data = dat[which(dat$year == 45),])
se.mreg45 <- sqrt(diag(vcovHC(mreg45)))


## Output to stargazer ##
stargazer(mreg1, mreg2, mreg42, mreg43, mreg44, mreg45,
          se = list(se.mreg1, se.mreg2, se.mreg42, se.mreg43, 
                    se.mreg44, se.mreg45),
          no.space = T, df = F,
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          dep.var.labels = c("Proportion Enlisted in County"),
          title = ("Impact of county presidential vote share on U.S. Army enlistment rates in World War 2. Robust standard errors listed in parentheses."),
          covariate.labels = c("Swing County", "New England", 
                               "Confederacy", "County GOP Vote Share",
                               "Pct. Urban", "Total Population (1000s)",
                               "Crop Value", "Pct. Male Eligible", "Pct. Black",
                               "Farms/Population",
                               "Mid Atlantic", "North East Central", 
                               "West North Central", "South Atlantic", 
                               "East South Central", "West South Central",
                               "Mountain", "Pacific",
                               "1942", "1943", "1944", "1945",
                               "Intercept"))


## Predicting Democratic Vote Share -- Table 3 Manuscript ##
swingm1 <- lm(demswing4044 ~ propelig + swingcounty
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop
              + as.factor(year), 
              data = dat[which(dat$year >= 41 & dat$year <= 43),])
se.swingm1 <- sqrt(diag(vcovHC(swingm1)))

swingm2 <- lm(demswing4044 ~ propelig * swingcounty
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop
              + as.factor(year), 
              data = dat[which(dat$year >= 41 & dat$year <= 43),])
se.swingm2 <- sqrt(diag(vcovHC(swingm2)))

swingm3 <- lm(demswing4044 ~ propelig + swingcounty
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop, 
              data = dat[which(dat$year == 43),])
se.swingm3 <- sqrt(diag(vcovHC(swingm3)))

# and output
stargazer(swingm1, swingm2, swingm3,
          se = list(se.swingm1, se.swingm2, se.swingm3),
          no.space = T, df = F,
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          dep.var.labels = c("Democratic Vote Share Swing (1944-1940"),
          title = c("Associations between enlistment rates and vote share swing towards the Democratic presidential candidate, 1944 election. Robust standard errors in parentheses."),
          covariate.labels = c("Enlistment Rate", "Swing County",
                               "Pct. Urban", "Population (1000s)",
                               "Avg. Crop Value", "Pct Male (21+)",
                               "Pct. Black", "Pct. Farms",
                               "1942", "1943",
                               "Enlistment Rate * Swing County",
                               "Intercept"))



## Table 1 -- Stepwise Regression
sr1 <- feols(propelig ~ swingcounty + yr42 + yr43 + yr44 + yr45,
                         data = dat.analysis,
                         cluster = "fips")
sr2 <- feols(propelig ~ swingcounty + electionyear +
               yr42 + yr43 + yr44 + yr45,
             data = dat.analysis,
             cluster = "fips")
sr3 <- feols(propelig ~ swingcounty + electionyear + GOPAvg +
               yr42 + yr43 + yr44 + yr45,
             data = dat.analysis,
             cluster = "fips")
sr4 <- feols(propelig ~ swingcounty + electionyear + GOPAvg + pcturban +
               yr42 + yr43 + yr44 + yr45,
             data = dat.analysis,
             cluster = "fips")
sr5 <- feols(propelig ~ swingcounty + electionyear + GOPAvg + pcturban +
             totpop.thou + 
             yr42 + yr43 + yr44 + yr45,
             data = dat.analysis,
             cluster = "fips")
sr6 <- feols(propelig ~ swingcounty + electionyear + GOPAvg + pcturban +
               totpop.thou + avgcropvalue +
               yr42 + yr43 + yr44 + yr45,
             data = dat.analysis,
             cluster = "fips")
sr7 <- feols(propelig ~ swingcounty + electionyear + GOPAvg + pcturban +
               totpop.thou + avgcropvalue + pctmale.elig2 + 
               yr42 + yr43 + yr44 + yr45,
             data = dat.analysis,
             cluster = "fips")
sr8 <- feols(propelig ~ swingcounty + electionyear + GOPAvg + pcturban +
               totpop.thou + avgcropvalue + pctmale.elig2 + pctblack +
               yr42 + yr43 + yr44 + yr45,
             data = dat.analysis,
             cluster = "fips")
sr9 <- feols(propelig ~ swingcounty + electionyear + GOPAvg + pcturban +
               totpop.thou + avgcropvalue + pctmale.elig2 + pctblack + farmprop + 
               yr42 + yr43 + yr44 + yr45,
             data = dat.analysis,
             cluster = "fips")

## output ##
etable(sr1, sr2, sr3, sr4, sr5, sr6, sr7, sr8, sr9,
       cluster = "fips", tex = TRUE,
       drop = "fips", digits = "r3", digits.stats = "r3",
       title = "Stepwise regression: Swing county coefficient when adding covariates.",
       label = "stepwisereg",
       fitstat = c("n", "ar2", "f"),
       dict = c(swingcounty = "Swing County (50-55)",
                GOPAvg = "GOP Vote Share",
                pcturban = "Pct Urban",
                totpop.thou = "Cty Pop (1000s)",
                avgcropvalue = "Avg Crop Value",
                pctmale.elig2 = "Pct Male Eligible",
                pctblack = "Pct Black",
                farmprop = "Prop. Farms",
                yr42 = "1942", yr43 = "1943", yr44 = "1944", yr45 = "1945",
                electionyear = "Election Year"))


## Table 2 -- Estimation with Gamma
dat.analysis$propelig.gamma <- dat.analysis$propelig + 0.00001

gamma1 <- glm(propelig.gamma ~ swingcounty
              + GOPAvg 
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig
              + pctblack + farmprop
              + as.factor(year), 
              data = dat.analysis, family = "Gamma"(link = "log"))
se.gamma1 <- sqrt(diag(vcovHC(gamma1)))
cl.gamma1 <- as.matrix(clusterMod(gamma1, dat.analysis$fips)[[1]])[,2]

gamma2 <- glm(propelig.gamma ~ swingcounty
              + electionyear + GOPAvg 
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig
              + pctblack + farmprop, 
              data = dat.analysis, family = "Gamma"(link = "log"))
se.gamma2 <- sqrt(diag(vcovHC(gamma2)))
cl.gamma2 <- as.matrix(clusterMod(gamma2, dat.analysis$fips)[[1]])[,2]

gamma3 <- glm(propelig.gamma ~ swingcounty
              * electionyear + GOPAvg 
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig
              + pctblack + farmprop, 
              data = dat.analysis, family = "Gamma"(link = "log"))
se.gamma3 <- sqrt(diag(vcovHC(gamma3)))
cl.gamma3 <- as.matrix(clusterMod(gamma3, dat.analysis$fips)[[1]])[,2]

stargazer(gamma1, gamma2, gamma3, 
          se = list(se.gamma1, se.gamma2, se.gamma3),
          no.space = T, df = F,
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          dep.var.labels = c("Proportion Enlisted in County"),
          title = ("Impact of congressional election results on U.S. Army enlistment rates in World War 2. Robust standard errors listed in parentheses."),
          covariate.labels = c("Swing County",
                               "Election Year", "County GOP Vote Share",
                               "Pct. Urban", "Total Population (1000s)",
                               "Crop Value", "Pct. Male Eligible", "Pct. Black",
                               "Farms/Population", "1942", "1943", "1944", "1945",
                               "Swing County x Election Year", "Intercept"))


## Table 3 -- State Fixed Effects
dat$electionyear <- ifelse(dat$year == 42 | dat$year == 44, 1, 0)

statefe1 <- lm(propelig ~ swingcounty
               + GOPAvg 
               + pcturban + pctmale.elig2
               + totpop.thou + avgcropvalue 
               + pctblack + farmprop
               + as.factor(year)
               + as.factor(state), data = dat)
se.statefe1 <- sqrt(diag(vcovHC(statefe1)))
cl.statefe1 <- as.matrix(clusterMod(statefe1, dat$fips)[[1]])[,2]

statefe2 <- lm(propelig ~ swingcounty + electionyear
               + GOPAvg 
               + pcturban + pctmale.elig2
               + totpop.thou + avgcropvalue 
               + pctblack + farmprop
               + as.factor(state), data = dat)
se.statefe2 <- sqrt(diag(vcovHC(statefe2)))
cl.statefe2 <- as.matrix(clusterMod(statefe2, dat$fips)[[1]])[,2]

statefe3 <- lm(propelig ~ swingcounty * electionyear
               + GOPAvg 
               + pcturban + pctmale.elig2
               + totpop.thou + avgcropvalue 
               + pctblack + farmprop
               + as.factor(state), data = dat)
se.statefe3 <- sqrt(diag(vcovHC(statefe3)))
cl.statefe3 <- as.matrix(clusterMod(statefe3, dat$fips)[[1]])[,2]

stargazer(statefe1, statefe2, statefe3,
          se = list(cl.statefe1, cl.statefe2, cl.statefe3),
          no.space = T, df = F, digits = 3,
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          omit = c("state"),
          title = ("Impact of congressional election results on U.S. Army enlistment rates in World War 2. Robust standard errors listed in parentheses."),
          covariate.labels = c("Swing County", "Election Year", "GOP Avg.", 
                               "Pct. Urban", "Pct. Male Eligible", "Total Pop. (1000s)",
                               "Avg. Crop Value", "Pct. Black", "Farms/Population",
                               "1941", "1942", "1943", "1944", "1945",
                               "Swing County x Election Year", "Intercept"))


## Table 4 -- Swing and Safe counties, swing and safe states
dat.analysis$state <- as.character(paste(dat.analysis$state, sep = ""))
dat.analysis$swingstate <- ifelse(dat.analysis$state == "massachusetts" |
                                    dat.analysis$state == "new hampshire" |
                                    dat.analysis$state == "michigan" | 
                                    dat.analysis$state == "wisconsin" | 
                                    dat.analysis$state == "illinois" | 
                                    dat.analysis$state == "new jersey" | 
                                    dat.analysis$state == "minnesota" | 
                                    dat.analysis$state == "new york" | 
                                    dat.analysis$state == "ohio" | 
                                    dat.analysis$state == "missouri", 1, 0)

# table of enlistment rates
tmp <- dat.analysis %>%
  group_by(swingstate, swingcounty) %>%
  summarise(mean.elig = mean(na.omit(propelig)))

# t-tests of enlistment rates
t1 <- t.test(propelig ~ swingstate, data = dat.analysis[which(dat.analysis$swingcounty == 1),])
t2 <- t.test(propelig ~ swingstate, data = dat.analysis[which(dat.analysis$swingcounty == 0),])
t3 <- t.test(propelig ~ swingcounty, data = dat.analysis[which(dat.analysis$swingstate == 1),])
t4 <- t.test(propelig ~ swingcounty, data = dat.analysis[which(dat.analysis$swingstate == 0),])

t.table <- matrix(NA, nrow = 5, ncol = 4)
t.table[2,2:4] <- c(round(t1$estimate, 4), round(t1$statistic, 4))
t.table[3,2:4] <- c(round(t2$estimate, 4), round(t2$statistic, 4))
t.table[4,2:4] <- c(round(t3$estimate, 4), round(t3$statistic, 4))
t.table[5,2:4] <- c(round(t4$estimate, 4), round(t4$statistic, 4))
t.table[1,] <- c("", "Est Group 1", "Est Group 2", "t-statistic")
t.table[,1] <- c("", "Swing Counties, Safe States/Swing States",
                 "Safe Counties, Safe States/Swing States",
                 "Swing States, Safe Counties/Swing Counties",
                 "Safe States, Safe Counties/Swing Counties")
stargazer(t.table)


## Table 5 -- Swing and Safe States
dat.swing <- dat.analysis %>% filter(swingstate == 1)
dat.safe <- dat.analysis %>% filter(swingstate == 0)
m.swingstate <- plm(propelig ~ swingcounty
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop
                    + as.factor(year),
                    data = dat.swing,
                    effect = "time", index = c("year"),
                    model = c("random"))
summary(m.swingstate) 
se.m.swingstate <- coeftest(m.swingstate, 
                            vcov=vcovHC(m.swingstate, type = "sss", cluster = "time"))[,2]
se.cl.swingstate <- coeftest(m.swingstate, 
                             vcov=vcovHC(m.swingstate, type = "sss", cluster = "group"))[,2]

# and safe state
m.safestate <- plm(propelig ~ swingcounty
                   + GOPAvg + pcturban + totpop.thou
                   + avgcropvalue + pctmale.elig2
                   + pctblack + farmprop
                   + as.factor(year),
                   data = dat.safe,
                   effect = "time", index = c("year"),
                   model = c("random"))
summary(m.safestate) 
se.m.safestate <- coeftest(m.safestate, 
                           vcov=vcovHC(m.safestate, type = "sss", cluster = "time"))[,2]
se.cl.safestate <- coeftest(m.safestate, 
                           vcov=vcovHC(m.safestate, type = "sss", cluster = "group"))[,2]

# and stargaze
stargazer(m.swingstate,
          m.safestate,
          se = list(se.m.swingstate,
                    se.m.safestate),
          no.space = T, df = F,
          omit = c("fips"),
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          dep.var.labels = c("Proportion Enlisted in County/Year"),
          title = ("Impact of county partisanship on U.S. Army enlistment rates in World War 2. Robust standard errors listed in parentheses."),
          covariate.labels = c("Swing County",
                               "County GOP Vote Share",
                               "Pct. Urban", "Total Population (1000s)",
                               "Crop Value", "Pct. Male Eligible", "Pct. Black",
                               "Farms/Population", 
                               "1942", "1943", "1944", "1945",
                               "Intercept")) 


## Table 6 --Congressional results -- at-large seats only ##

# t-tests
dat.analysis <- dat.analysis %>%
  mutate(cong_atlarge = ifelse(state == "arizona" | state == "delaware" | state == "nevada" | state == "new mexico" | state == "north dakota" | state == "vermont" | state == "wyoming" | state == "connecticut" | state == "florida " | state == "illinois" | state == "new york" | state == "ohio" | state == "pennsylvania", 1, 0))

tc1 <- t.test(propelig ~ swingcongavg, data = dat.analysis[which(dat.analysis$cong_atlarge == 1),])
tc2 <- t.test(propelig ~ swingcongavg, data = dat.analysis[which(dat.analysis$cong_atlarge == 0),])

# Output
tcs <- matrix(NA, ncol = 4, nrow = 3)
tcs[1,] <- c("", "Est Group 1", "Est Group 2", "t-statistic")
tcs[2,] <- c("At-Large CDs, Safe Counties/Swing Counties", "0.911", "0.840", "1.90")
tcs[3,] <- c("Not At-Large CDs, Safe Counties/Swing Counties", "0.929", "0.858", "3.30")
xtable(tcs)

# Regression results
cong_atlarge <- dat.analysis %>%
  filter(state == "arizona" | state == "delaware" | state == "nevada" | state == "new mexico" | state == "north dakota" | state == "vermont" | state == "wyoming" | state == "connecticut" | state == "florida " | state == "illinois" | state == "new york" | state == "ohio" | state == "pennsylvania")

cong_notatlarge <- dat.analysis %>%
  filter(state != "arizona"  | state != "delaware" | state != "nevada" | state != "new mexico" | state != "north dakota" | state != "vermont" | state != "wyoming" | state != "connecticut" | state != "florida " | state != "illinois" | state != "new york" | state != "ohio" | state != "pennsylvania")

mc_atlarge <- feols(propelig ~ swingcongavg
                    + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop 
                    + yr42 + yr43 + yr44 + yr45,
                    data = cong_atlarge,
                    cluster = "fips")

mc_notatlarge <- feols(propelig ~ swingcongavg
                       + pcturban + totpop.thou
                       + avgcropvalue + pctmale.elig2
                       + pctblack + farmprop 
                       + yr42 + yr43 + yr44 + yr45,
                       data = cong_notatlarge,
                       cluster = "fips")

# Output
etable(mc_atlarge, mc_notatlarge,
       cluster = "fips", tex = TRUE,
       drop = "statefips", digits = "r3", digits.stats = "r3",
       title = "Congressional analyses, statewide districts only.",
       label = "congstatewide",
       fitstat = c("n", "ar2", "f"))


## Table 7 -- Volunteers and draftees
# volunteers
m1vol <- plm(propvol ~ swingcounty
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop
             + as.factor(year),
             data = dat.analysis,
             effect = "time", index = c("year"),
             model = c("random"))
summary(m1vol) 
se.m1vol <- coeftest(m1vol, vcov=vcovHC(m1vol, type = "sss", cluster = "time"))[,2]
cl.m1vol <- coeftest(m1vol, vcov=vcovHC(m1vol, type = "sss", cluster = "group"))[,2]

m2vol <- lm(propvol ~ swingcounty + electionyear
            + GOPAvg + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat.analysis)
se.m2vol <- sqrt(diag(vcovHC(m2vol)))
cl.m2vol <- as.matrix(clusterMod(m2vol, dat.analysis$fips)[[1]])[,2]

m3vol <- lm(propvol ~ swingcounty * electionyear
            + GOPAvg + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat.analysis)
se.m3vol <- sqrt(diag(vcovHC(m3vol)))
cl.m3vol <- as.matrix(clusterMod(m3vol, dat.analysis$fips)[[1]])[,2]

# draftees
m1dra <- plm(propdraft ~ swingcounty
             + GOPAvg + pcturban + totpop.thou
             + avgcropvalue + pctmale.elig2
             + pctblack + farmprop
             + as.factor(year),
             data = dat[which(dat$year >= 41),],
             effect = "time", index = c("year"),
             model = c("random"))
se.m1dra <- coeftest(m1dra, vcov=vcovHC(m1dra, type = "sss", cluster = "time"))[,2]
cl.m1dra <- coeftest(m1dra, vcov=vcovHC(m1dra, type = "sss", cluster = "group"))[,2]

m2dra <- lm(propdraft ~ swingcounty + electionyear
            + GOPAvg + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat.analysis)
se.m2dra <- sqrt(diag(vcovHC(m2dra)))
cl.m2dra <- as.matrix(clusterMod(m2dra, dat.analysis$fips)[[1]])[,2]

m3dra <- lm(propdraft ~ swingcounty * electionyear
            + GOPAvg + pcturban + totpop.thou
            + avgcropvalue + pctmale.elig2
            + pctblack + farmprop,
            data = dat.analysis)
se.m3dra <- sqrt(diag(vcovHC(m3dra)))
cl.m3dra <- as.matrix(clusterMod(m3dra, dat.analysis$fips)[[1]])[,2]

## Output to stargazer ##
stargazer(m1vol, m2vol, m3vol,
          m1dra, m2dra, m3dra,
          se = list(se.m1vol, se.m2vol, se.m3vol,
                    se.m1dra, se.m2dra, se.m3dra),
          no.space = T, df = F,
          omit = c("state", "fips"),
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          dep.var.labels = c("Volunteers", "Draftees"),
          title = ("Impact of county partisanship on U.S. Army enlistment rates in World War 2. Robust standard errors listed in parentheses."),
          covariate.labels = c("Swing County",
                               "Election Year", "County GOP Vote Share",
                               "Pct. Urban", "Total Population (1000s)",
                               "Crop Value", "Pct. Male Eligible", "Pct. Black",
                               "Farms/Population", "1942", "1943", "1944", "1945",
                               "Swing County x Election Year", "Intercept"))


## Table 8 -- Matching as pre-processing
t.test.out <- matrix(NA, nrow = 400, ncol = 6)

for(i in 1:400){
  
  tmp <- as.data.frame(cbind(dat.census$propelig, 
                             dat.census$swingcounty,
                             dat.census[,i+85]))
  if(is.numeric(tmp$V3) == F) next
  if(min(na.omit(tmp$V3)) == max(na.omit(tmp$v3))) next
  
  test <- t.test(tmp$V3 ~ tmp$V2)
  
  t.test.out[i,2] <- as.numeric(paste(test$estimate[1], sep = ""))
  t.test.out[i,3] <- as.numeric(paste(test$estimate[2], sep = ""))
  t.test.out[i,4] <- as.numeric(paste(test$statistic, sep = ""))
  t.test.out[i,5] <- as.numeric(paste(test$p.value, sep = ""))
  t.test.out[i,6] <- as.numeric(paste(test$parameter, sep = ""))
  
  cat("\r", i, "of", 400)
}

t.test.out[,1] <- colnames(dat.census)[86:485]
t.test.out <- as.data.frame(t.test.out)
colnames(t.test.out) <- c("test.var", "group1", "group2", "tstat", "pvalue", "df")
t.test.out$pvalue <- as.numeric(paste(t.test.out$pvalue, sep = ""))

## Find where pvalue <= 0.05 ##
different.counties <- t.test.out[which(t.test.out$pvalue <= 0.05),]

## create dataset of theoretically-derived vars from different.counties ##
prematch.dat <- dat.census[,c("propelig", "swingcounty", "electionyear",
                              "GOPAvg", "pcturban", "totpop.thou",
                              "avgcropvalue", "pctmale.elig2", "pctblack",
                              "farmprop", "year", "tsch1820", "tsch2124",
                              "m25edu0", "m25edu14", "m25col13", "m25col4",
                              "rurfarm", "pctm14lf", "mfarm", "mlabor",
                              "wholwage.y", "servwage.y", "pcturb40",
                              "medperre", "pctradio", "pctrefri", "avrent",
                              "meschm25", "unempnem", "density", "fips")]

# Define vars
Y = prematch.dat$propelig
Tr = prematch.dat$swingcounty
X = prematch.dat[,c(3:32)]

t <- na.omit(cbind(Y, Tr, X))

## Estimate balance on Tr ##
tr.predict <- glm(Tr ~ ., data = prematch.dat[, -c(1, 32)])
summary(tr.predict)

Z <- cbind(t[,c(4:6, 8:11, 16:24, 26, 27, 30)])

## Match ##
match <- Match(Y=t$Y, Tr=t$Tr, X=t[,c(3:31)], 
               Z=Z, BiasAdjust=T, estimand="ATT", 
               M=1, replace=T)
summary(match)

# check for balance 
postmatch.dat <- as.data.frame(rbind(t[match$index.treated,],
                                     t[match$index.control,]))
tr.postdict <- glm(Tr ~ ., data = postmatch.dat[, -c(1)])
summary(tr.postdict)

## Regression, m1 from manuscript
library(lmtest)
match.fe1 <- lm(Y ~ Tr + as.factor(year) + as.factor(fips),
          data = postmatch.dat)
se.match.fe1 <- as.matrix(clusterMod(match.fe1, postmatch.dat$fips)[[1]])[,2]

match.fe2 <- lm(Y ~ Tr * as.factor(year) + as.factor(fips), data = postmatch.dat)
se.match.fe2 <- as.matrix(clusterMod(match.fe2, postmatch.dat$fips)[[1]])[,2]

match.m1 <- plm(Y ~ Tr
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop
                + as.factor(year),
                data = postmatch.dat,
                effect = "time", index = c("year"),
                model = c("random"))
se.match.m1 <- coeftest(match.m1, 
                        vcov=vcovHC(match.m1, 
                                        type = "sss", 
                                        cluster = "time"))[,2]

match.m2 <- plm(Y ~ Tr + electionyear
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop,
                data = postmatch.dat,
                effect = "time", index = c("year"),
                model = c("random"))
se.match.m2 <- coeftest(match.m2, 
                        vcov=vcovHC(match.m2, 
                                        type = "sss", 
                                        cluster = "time"))[,2]

match.m3 <- plm(Y ~ Tr * electionyear
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop,
                data = postmatch.dat,
                effect = "time", index = c("year"),
                model = c("random"))
se.match.m3 <- coeftest(match.m3, 
                        vcov=vcovHC(match.m3, 
                                        type = "sss", 
                                        cluster = "time"))[,2]

## Output to stargazer ##
stargazer(match.fe1, match.fe2,
          match.m1, match.m2, match.m3,
          se = list(se.match.fe1, se.match.fe2,
                    se.match.m1, se.match.m2, se.match.m3),
          no.space = T, df = F, omit = c("fips"),
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          dep.var.labels = c("Proportion Enlisted in County/Year"),
          title = ("Impact of county partisanship on U.S. Army enlistment rates in World War 2. Robust standard errors listed in parentheses."),
          covariate.labels = c("Swing County",
                               "Election Year", "County GOP Vote Share",
                               "Pct. Urban", "Total Population (1000s)",
                               "Crop Value", "Pct. Male Eligible", "Pct. Black",
                               "Farms/Population", "1942", "1943", "1944", "1945",
                               "Swing County x 1942", "Swing County x 1943", 
                               "Swing County x 1944", "Swing County x 1945", 
                               "Swing County x Election Year", "Intercept"))

## Create covariate balance plot ##
library(cobalt); library(MatchIt)
balance.table <- bal.tab(match, treat = t$Tr, covs = t[,c(3:31)],
                         estimand = "ATT")
covs <- t[,c(3:31)]

# and plot
m.out <- matchit(f.build("Tr", covs), data = t,
                 method = "nearest", replace = T)
love.out <- love.plot(m.out, thresholds = c(m = 0.1, ks = 0.1),
                      stats = c("m", "ks.statistics"), binary = "std", abs = T,
                      var.order = "unadjusted", sample.names = c("Unmatched", "Matched"),
                      position = "top", shapes = c("circle filled", "diamond filled"), 
                      colors = c("mediumorchid4", "goldenrod2"))
png("loveplot.png", res = 100)
print(love.out)
dev.off() # not working; manually exported as .png


# original/baseline models, unit fixed effects
utfe_all <- feols(propelig ~ swingcounty
                  + as.factor(year) + as.factor(fips),
                  data = dat.analysis,
                  cluster = "fips")
utfe_42 <- feols(propelig ~ swingcounty
                 + as.factor(fips), 
                 data = dat42,
                 se = "hetero")
utfe_43 <- feols(propelig ~ swingcounty
                 + as.factor(fips), 
                 data = dat43,
                 se = "hetero")
utfe_44 <- feols(propelig ~ swingcounty
                 + as.factor(fips), 
                 data = dat44,
                 se = "hetero")
utfe_45 <- feols(propelig ~ swingcounty
                 + as.factor(fips), 
                 data = dat45,
                 se = "hetero")
utfe_am1 <- feols(propelig ~ splitswing
                  + as.factor(year) + as.factor(fips), 
                  data = dat.analysis,
                  cluster = "fips")
utfe_am2 <- feols(propelig ~ pickupswing
                  + as.factor(year) + as.factor(fips), 
                  data = dat.analysis,
                  cluster = "fips")

etable(utfe_all, utfe_42, utfe_43, utfe_44, utfe_45, utfe_am1, utfe_am2,
       cluster = "fips", tex = TRUE,
       drop = c("fips", "year"), digits = "r3", digits.stats = "r3",
       title = "Impact of county partisanship on U.S. Army enlistment rates in World War 2. Unit-clustered standard errors listed in parentheses.",
       label = "reg1",
       fitstat = c("n", "ar2", "f"),
       dict = c(swingcounty = "Swing County (50-55)",
                GOPAvg = "GOP Vote Share",
                pcturban = "Pct Urban",
                totpop.thou = "Cty Pop (1000s)",
                avgcropvalue = "Avg Crop Value",
                pctmale.elig2 = "Pct Male Eligible",
                pctblack = "Pct Black",
                farmprop = "Prop. Farms",
                yr42 = "1942", yr43 = "1943", yr44 = "1944", yr45 = "1945",
                electionyear = "Election Year",
                splitswing = "Swing County (47.5-52.5)",
                pickupswing = "Swing County (45-50)"))


## Table 13 -- Checking Linear Interpolation
pop <- dat %>%
  group_by(state, year) %>%
  summarise(statepop = sum(totpop.interpolate)) %>%
  mutate(lagpop = dplyr::lag(statepop, n = 1, default = NA))
pop <- left_join(pop, vitals, by = c("state", "year"))

# Regress births and deaths on state population
pop.m1 <- lm(statepop ~ births + deaths, data = pop)
pop.m2 <- lm(statepop ~ births + deaths + lagpop, data = pop)
stargazer(pop.m1, pop.m2)


## Table 14: Volunteers Pre-war
prewar.dat <- dat %>% filter(year <= 41)

ws0 <- cor(dat.analysis$propvol, dat.analysis$propdraft,
           use = "pairwise.complete.obs")
wss <- cor(prewar.dat$propvol, prewar.dat$swingcounty,
           use = "pairwise.complete.obs")

ws1 <- glm(swingcounty ~ propvol,
           data = prewar.dat,
           family = binomial(link = "logit"))
se.ws1 <- as.matrix(clusterMod(ws1, prewar.dat$fips)[[1]])[,2]

ws2 <- glm(swingcounty ~ propvol
           + GOPAvg + pcturban + totpop.thou
           + avgcropvalue + pctmale.elig2
           + pctblack + farmprop,
           data = prewar.dat,
           family = binomial(link = "logit"))
se.ws2 <- as.matrix(clusterMod(ws2, prewar.dat$fips)[[1]])[,2]

dat.4142 <- dat %>% filter(year == 41 | year == 42) %>%
  mutate(propvol.lag = dplyr::lag(propvol, n = 1, default = NA),
         propelig.lag = dplyr::lag(propelig, n = 1, default = NA),
         propdraft.lag = dplyr::lag(propdraft, n = 1, default = NA))

ws3 <- lm(propelig ~ swingcounty + propvol.lag + propdraft.lag
          + GOPAvg + pcturban + totpop.thou
          + avgcropvalue + pctmale.elig2
          + pctblack + farmprop,
          data = dat.4142)
se.ws3 <- sqrt(diag(vcovHC(ws3)))

# output
stargazer(ws1, ws2, ws3,
          se = list(se.ws1, se.ws2, se.ws3),
          no.space = T, df = F,
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          dep.var.labels = c("Swing County Status"),
          title = ("Association between volunteer rates, US army, 1940 and 1941, and New Deal-era swing county status."),
          covariate.labels = c("Pct Volunteers", "Swing County",
                               "Pct Volunteers_{t-1}", "Pct Enlistments_{t-1}",
                               "GOP Avg", "Pct Urban",
                               "Total Pop. (1000s)", "Avg. Crop Value",
                               "Pct. Male Eligible", "Pct. Black", "Farms/Pop",
                               "Intercept"))


##### Manuscript extension -- Predicting 1944 presidential election #####
intent1 <- lm(demswing4044 ~ propelig + swingcounty
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop
              + as.factor(year), 
              data = dat.analysis)
se.intent1 <- sqrt(diag(vcovHC(intent1)))

intent2 <- lm(demswing4044 ~ propelig * swingcounty
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop
              + as.factor(year), 
              data = dat.analysis)
se.intent2 <- sqrt(diag(vcovHC(intent2)))

intent3 <- lm(demswing4044 ~ propelig + swingcounty
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop, 
              data = dat.analysis)
se.intent3 <- sqrt(diag(vcovHC(intent3)))

# and output
stargazer(intent1, intent2, intent3,
          se = list(se.intent1, se.intent2, se.intent3),
          no.space = T, df = F,
          keep.stat = c("n", "adj.rsq", "f", "ll", "aic", "theta"),
          dep.var.labels = c("Democratic Vote Share Swing (1944-1940"),
          title = c("Associations between enlistment rates and vote share swing towards the Democratic presidential candidate, 1944 election. Robust standard errors in parentheses."),
          covariate.labels = c("Enlistment Rate", "Swing County",
                               "Pct. Urban", "Population (1000s)",
                               "Avg. Crop Value", "Pct Male (21+)",
                               "Pct. Black", "Pct. Farms",
                               "1942", "1943", "1944", "1945",
                               "Enlistment Rate * Swing County",
                               "Intercept"))


## Congressional models, all dependent variables
# original/baseline models
cong1 <- feols(propelig ~ swingcongavg
               + GOPAvg + pcturban + totpop.thou
               + avgcropvalue + pctmale.elig2
               + pctblack + farmprop
               + yr42 + yr43 + yr44 + yr45,
               data = dat.analysis,
               cluster = "fips")
cong2 <- feols(propelig ~ swingcong36
               + GOPAvg + pcturban + totpop.thou
               + avgcropvalue + pctmale.elig2
               + pctblack + farmprop,
               data = dat42,
               cluster = "fips")
cong3 <- feols(propelig ~ splitcongswing36
               + GOPAvg + pcturban + totpop.thou
               + avgcropvalue + pctmale.elig2
               + pctblack + farmprop,
               data = dat42,
               cluster = "fips")
cong4 <- feols(propelig ~ pickupcongswing36
               + GOPAvg + pcturban + totpop.thou
               + avgcropvalue + pctmale.elig2
               + pctblack + farmprop, 
               data = dat42,
               se = "hetero")
cong5 <- feols(propelig ~ swingcong38
               + GOPAvg + pcturban + totpop.thou
               + avgcropvalue + pctmale.elig2
               + pctblack + farmprop, 
               data = dat42,
               se = "hetero")
cong6 <- feols(propelig ~ splitcongswing38
               + GOPAvg + pcturban + totpop.thou
               + avgcropvalue + pctmale.elig2
               + pctblack + farmprop, 
               data = dat42,
               se = "hetero")
cong7 <- feols(propelig ~ pickupcongswing38
               + GOPAvg + pcturban + totpop.thou
               + avgcropvalue + pctmale.elig2
               + pctblack + farmprop, 
               data = dat42,
               se = "hetero")
cong8 <- feols(propelig ~ swingcong40
               + GOPAvg + pcturban + totpop.thou
               + avgcropvalue + pctmale.elig2
               + pctblack + farmprop, 
               data = dat42,
               cluster = "fips")
cong9 <- feols(propelig ~ splitcongswing40
               + GOPAvg + pcturban + totpop.thou
               + avgcropvalue + pctmale.elig2
               + pctblack + farmprop, 
               data = dat42,
               cluster = "fips")
cong10 <- feols(propelig ~ pickupcongswing40
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong11 <- feols(propelig ~ swingcong42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                cluster = "fips")
cong12 <- feols(propelig ~ splitcongswing42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                cluster = "fips")
cong13 <- feols(propelig ~ pickupcongswing42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                se = "hetero")

# propdraft variant
cong14 <- feols(propdraft ~ swingcongavg
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop
                + yr42 + yr43 + yr44 + yr45,
                data = dat.analysis,
                cluster = "fips")
cong15 <- feols(propdraft ~ swingcong36
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop,
                data = dat42,
                cluster = "fips")
cong16 <- feols(propdraft ~ splitcongswing36
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop,
                data = dat42,
                cluster = "fips")
cong17 <- feols(propdraft ~ pickupcongswing36
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong18 <- feols(propdraft ~ swingcong38
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong19 <- feols(propdraft ~ splitcongswing38
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong20 <- feols(propdraft ~ pickupcongswing38
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong21 <- feols(propdraft ~ swingcong40
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                cluster = "fips")
cong22 <- feols(propdraft ~ splitcongswing40
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                cluster = "fips")
cong23 <- feols(propdraft ~ pickupcongswing40
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong24 <- feols(propdraft ~ swingcong42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                cluster = "fips")
cong25 <- feols(propdraft ~ splitcongswing42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                cluster = "fips")
cong26 <- feols(propdraft ~ pickupcongswing42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                se = "hetero")

# male-eligible interpolated DV variant
cong27 <- feols(male1554.elig ~ swingcongavg
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop
                + yr42 + yr43 + yr44 + yr45,
                data = dat.analysis,
                cluster = "fips")
cong28 <- feols(male1554.elig ~ swingcong36
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop,
                data = dat42,
                cluster = "fips")
cong29 <- feols(male1554.elig ~ splitcongswing36
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop,
                data = dat42,
                cluster = "fips")
cong30 <- feols(male1554.elig ~ pickupcongswing36
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong31 <- feols(male1554.elig ~ swingcong38
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong32 <- feols(male1554.elig ~ splitcongswing38
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong33 <- feols(male1554.elig ~ pickupcongswing38
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong34 <- feols(male1554.elig ~ swingcong40
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                cluster = "fips")
cong35 <- feols(male1554.elig ~ splitcongswing40
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                cluster = "fips")
cong36 <- feols(male1554.elig ~ pickupcongswing40
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong37 <- feols(male1554.elig ~ swingcong42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                cluster = "fips")
cong38 <- feols(male1554.elig ~ splitcongswing42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                cluster = "fips")
cong39 <- feols(male1554.elig ~ pickupcongswing42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                se = "hetero")

# draft_melig variant
cong40 <- feols(draft_melig ~ swingcongavg
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop
                + yr42 + yr43 + yr44 + yr45,
                data = dat.analysis,
                cluster = "fips")
cong41 <- feols(draft_melig ~ swingcong36
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop,
                data = dat42,
                cluster = "fips")
cong42 <- feols(draft_melig ~ splitcongswing36
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop,
                data = dat42,
                cluster = "fips")
cong43 <- feols(draft_melig ~ pickupcongswing36
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong44 <- feols(draft_melig ~ swingcong38
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong45 <- feols(draft_melig ~ splitcongswing38
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong46 <- feols(draft_melig ~ pickupcongswing38
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong47 <- feols(draft_melig ~ swingcong40
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                cluster = "fips")
cong48 <- feols(draft_melig ~ splitcongswing40
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                cluster = "fips")
cong49 <- feols(draft_melig ~ pickupcongswing40
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat42,
                se = "hetero")
cong50 <- feols(draft_melig ~ swingcong42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                cluster = "fips")
cong51 <- feols(draft_melig ~ splitcongswing42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                cluster = "fips")
cong52 <- feols(draft_melig ~ pickupcongswing42
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop, 
                data = dat43,
                se = "hetero")

# and tidy
tidycong <- rbind(tidy(cong1) %>%
                    filter(term == "swingcongavg"),
                  tidy(cong2) %>%
                    filter(term == "swingcong36"),
                  tidy(cong3) %>%
                    filter(term == "splitcongswing36"),
                  tidy(cong4) %>%
                    filter(term == "pickupcongswing36"),
                  tidy(cong5) %>%
                    filter(term == "swingcong38"),
                  tidy(cong6) %>%
                    filter(term == "splitcongswing38"),
                  tidy(cong7) %>%
                    filter(term == "pickupcongswing38"),
                  tidy(cong8) %>%
                    filter(term == "swingcong40"),
                  tidy(cong9) %>%
                    filter(term == "splitcongswing40"),
                  tidy(cong10) %>%
                    filter(term == "pickupcongswing40"),
                  tidy(cong11) %>%
                    filter(term == "swingcong42"),
                  tidy(cong12) %>%
                    filter(term == "splitcongswing42"),
                  tidy(cong13) %>%
                    filter(term == "pickupcongswing42"),
                  tidy(cong14) %>%
                    filter(term == "swingcongavg"),
                  tidy(cong15) %>%
                    filter(term == "swingcong36"),
                  tidy(cong16) %>%
                    filter(term == "splitcongswing36"),
                  tidy(cong17) %>%
                    filter(term == "pickupcongswing36"),
                  tidy(cong18) %>%
                    filter(term == "swingcong38"),
                  tidy(cong19) %>%
                    filter(term == "splitcongswing38"),
                  tidy(cong20) %>%
                    filter(term == "pickupcongswing38"),
                  tidy(cong21) %>%
                    filter(term == "swingcong40"),
                  tidy(cong22) %>%
                    filter(term == "splitcongswing40"),
                  tidy(cong23) %>%
                    filter(term == "pickupcongswing40"),
                  tidy(cong24) %>%
                    filter(term == "swingcong42"),
                  tidy(cong25) %>%
                    filter(term == "splitcongswing42"),
                  tidy(cong26) %>%
                    filter(term == "pickupcongswing42"),
                  tidy(cong27) %>%
                    filter(term == "swingcongavg"),
                  tidy(cong28) %>%
                    filter(term == "swingcong36"),
                  tidy(cong29) %>%
                    filter(term == "splitcongswing36"),
                  tidy(cong30) %>%
                    filter(term == "pickupcongswing36"),
                  tidy(cong31) %>%
                    filter(term == "swingcong38"),
                  tidy(cong32) %>%
                    filter(term == "splitcongswing38"),
                  tidy(cong33) %>%
                    filter(term == "pickupcongswing38"),
                  tidy(cong34) %>%
                    filter(term == "swingcong40"),
                  tidy(cong35) %>%
                    filter(term == "splitcongswing40"),
                  tidy(cong36) %>%
                    filter(term == "pickupcongswing40"),
                  tidy(cong37) %>%
                    filter(term == "swingcong42"),
                  tidy(cong38) %>%
                    filter(term == "splitcongswing42"),
                  tidy(cong39) %>%
                    filter(term == "pickupcongswing42"),
                  tidy(cong40) %>%
                    filter(term == "swingcongavg"),
                  tidy(cong41) %>%
                    filter(term == "swingcong36"),
                  tidy(cong42) %>%
                    filter(term == "splitcongswing36"),
                  tidy(cong43) %>%
                    filter(term == "pickupcongswing36"),
                  tidy(cong44) %>%
                    filter(term == "swingcong38"),
                  tidy(cong45) %>%
                    filter(term == "splitcongswing38"),
                  tidy(cong46) %>%
                    filter(term == "pickupcongswing38"),
                  tidy(cong47) %>%
                    filter(term == "swingcong40"),
                  tidy(cong48) %>%
                    filter(term == "splitcongswing40"),
                  tidy(cong49) %>%
                    filter(term == "pickupcongswing40"),
                  tidy(cong50) %>%
                    filter(term == "swingcong42"),
                  tidy(cong51) %>%
                    filter(term == "splitcongswing42"),
                  tidy(cong52) %>%
                    filter(term == "pickupcongswing42")
) %>%
  mutate(model = c(rep("enlist_pop", 13), rep("draft_pop", 13),
         rep("enlist_melig", 13), rep("draft_melig", 13)))

# coefplot output
library(dotwhisker)
cong_allmodels <- dwplot(tidycong, dodge_size = 0.75, ci = 0.95,
                         dot_args = list(aes(shape = model)),
                         whisker_args = list(aes(linetype = model)),
                         vline = geom_vline(xintercept = 0,
                                            colour = "grey60",
                                            linetype = 2)) +
  theme_bw() + xlab("Coefficient Estimate") + ylab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  ggtitle("") +
  theme(plot.title = element_text(face = "bold"),
        legend.title = element_blank(),
        legend.position = "bottom") + 
  guides(colour = guide_legend(override.aes = list(shape = c(0, 1, 2, 3), 
                                                   linetype = c(0, 1, 2, 3)))) + 
  scale_shape_manual(name = "Model",
                     values = c(0, 1, 2, 3),
                     guide = "none") +
  scale_linetype_manual(name = "Model",
                        values = c(0, 1, 2, 3),
                        guide = "none") + 
  scale_colour_grey(start = .2, end = .25,
                    guide = "legend")
png("cong_allmodels.png", res = 80)
print(cong_allmodels)
dev.off()



## Punitive Conscription
punish1 <- feols(propelig ~ demswing4044 
          + pcturban + totpop.thou
          + avgcropvalue + pctmale.elig2
          + pctblack + farmprop, data = dat45,
          cluster = "fips")

punish2 <- feols(draft_melig ~ demswing4044 
           + pcturban + totpop.thou
           + avgcropvalue + pctmale.elig2
           + pctblack + farmprop, data = dat45,
           cluster = "fips")

dat43$congswing <- dat43$demcong1942 - dat43$demcong1940

punish3 <- feols(propelig ~ congswing 
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop, dat = dat43,
              cluster = "fips")

punish4 <- feols(draft_melig ~ congswing 
              + pcturban + totpop.thou
              + avgcropvalue + pctmale.elig2
              + pctblack + farmprop, dat = dat43,
              cluster = "fips")

etable(punish1, punish2, punish3, punish4,
       cluster = "fips", tex = TRUE,
       drop = c("fips", "year"), digits = "r3", digits.stats = "r3",
       title = "Impact of county partisanship on U.S. Army enlistment rates in World War 2. Unit-clustered standard errors listed in parentheses.",
       label = "punitive",
       fitstat = c("n", "ar2", "f"),
       dict = c(demswing4044 = "Swing to Democrats (Pres.)",
                pcturban = "Pct Urban",
                totpop.thou = "Cty Pop (1000s)",
                avgcropvalue = "Avg Crop Value",
                pctmale.elig2 = "Pct Male Eligible",
                pctblack = "Pct Black",
                farmprop = "Prop. Farms",
                congswing = "Swing to Democrats (Cong.)"))


## WI/MI/IL Only
midwest <- dat.analysis %>%
  filter(state == "illinois" | state == "michigan" | state == "wisconsin")
midwest2 <- dat42 %>%
  filter(state == "illinois" | state == "michigan" | state == "wisconsin")

midwestm1 <-  feols(propelig ~ swingcounty 
                + GOPAvg + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop
                + yr42 + yr43 + yr44 + yr45, data = midwest,
                cluster = "fips")
midwestm2 <-  feols(propelig ~ swingcounty 
                    + GOPAvg + pcturban + totpop.thou
                    + avgcropvalue + pctmale.elig2
                    + pctblack + farmprop, data = midwest2,
                    cluster = "fips")
etable(midwestm1, midwestm2,
       cluster = "fips", tex = TRUE,
       drop = c("fips", "year"), digits = "r3", digits.stats = "r3",
       title = "Impact of county partisanship on U.S. Army enlistment rates in World War 2. Illinois, Michigan, and Wisconsin only. Unit-clustered standard errors listed in parentheses.",
       label = "midwestonly",
       fitstat = c("n", "ar2", "f"),
       dict = c(swingcounty = "Swing County",
                DemAvg = "Democratic Vote Share",
                GOPAvg = "GOP Vote Share",
                pcturban = "Pct Urban",
                totpop.thou = "Cty Pop (1000s)",
                avgcropvalue = "Avg Crop Value",
                pctmale.elig2 = "Pct Male Eligible",
                pctblack = "Pct Black",
                farmprop = "Prop. Farms",
                yr42 = "1942", yr43 = "1943", yr44 = "1944", yr45 = "1945"))


## DemAvg vs GOPAvg
demavg <- feols(propelig ~ swingcounty + DemAvg
                 + pcturban + totpop.thou
                 + avgcropvalue + pctmale.elig2
                 + pctblack + farmprop 
                 + yr42 + yr43 + yr44 + yr45, data = dat.analysis,
                 cluster = "fips")
gopavg <- feols(propelig ~ swingcounty + GOPAvg
                + pcturban + totpop.thou
                + avgcropvalue + pctmale.elig2
                + pctblack + farmprop 
                + yr42 + yr43 + yr44 + yr45, data = dat.analysis,
                cluster = "fips")
etable(demavg, gopavg,
       cluster = "fips", tex = TRUE,
       drop = c("fips", "year"), digits = "r3", digits.stats = "r3",
       title = "Impact of county partisanship on U.S. Army enlistment rates in World War 2. Unit-clustered standard errors listed in parentheses.",
       label = "demavg",
       fitstat = c("n", "ar2", "f"),
       dict = c(swingcounty = "Swing County",
                DemAvg = "Democratic Vote Share",
                GOPAvg = "GOP Vote Share",
                pcturban = "Pct Urban",
                totpop.thou = "Cty Pop (1000s)",
                avgcropvalue = "Avg Crop Value",
                pctmale.elig2 = "Pct Male Eligible",
                pctblack = "Pct Black",
                farmprop = "Prop. Farms",
                yr42 = "1942", yr43 = "1943", yr44 = "1944", yr45 = "1945"))
